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Section  1 
INTRODUCTION 


The  goal  of  the  theoretical  investigation  described  here  is  to  obtain  a  clearer 
understanding  of  the  fluid  dynamics  of  liquid-filled  shells  for  those  cases  where  the 
fluid  only  partially  fills  its  cavity.  Problems  can  arise  during  the  flight  of  such 
projectiles  when  there  is  a  strong  coupling  between  the  perturbed  motion  of  the  rapidly 
spinning  shell  casing  and  the  attendant  wave  motion  excited  in  the  fluid.  In  addition, 
unlike  a  solid  filler  which  always  rotates  with  the  same  angular  velocity  as  the  casing, 
the  liquid  filler  will  rotate  at  a  reduced  angular  velocity  over  much  of  the  trajectory. 
The  result  is  a  net  transfer  of  angular  momentum  from  casing  to  fluid,  which  in  itself 
can  be  destabilizing.1  Considerations  of  internal  fluid  motion  also  play  an  important 
role  in  spin-stabilized  liquid  propellant  rockets.2 

The  model  problem  most  often  studied  by  investigators  is  the  ax isym metric  spin- 
up  of  fluid  in  a  right  circular  cylinder.  We  seek  to  represent  the  evolution  with  time 
of  the  fluid  motion,  particularly  the  azimuthal  velocity,  as  it  asymptotically  approaches 
solid-body  rotation.  Such  a  calculation  is  necessary  as  input  for  the  prediction  of  the 
spin  decay  of  liquid-filled  projectiles.3’4  It  can  also  serve  as  the  base  flow  about  which 
to  carry  out  a  perturbation  analysis  of  the  eigenfrequencies  of  the  fluid.5’6 7  The  latter 
are  of  critical  importance  to  the  question  of  the  stability  of  the  fluid-shell  system.1 

Most  previous  investigations  of  the  problem  have  assumed  the  cylinder  to  be 
completely  filled.  The  earliest  work  appears  to  be  that  of  McLeod, ^  who  assumed 


1.  Engineering  Design  Handbook.  Liquid-Filled  Projectile  Design,  AMC  Pamphlet 
NoT  706-165,  05.  Army  Materiel  Command,  Washington,  D.C.,  April  1969. 

2.  Winch,  D.M.,  An  Investigation  of  the  Liquid  Level  at  the  Wall  of  a  Spinning 
Tank,  NASA  TN  D-1536,  August  1962. 

3.  Kitchens,  C.W.,  Jr.,  Gerber,  N.  and  Sedney,  R.,  "Spin  Decay  of  Liquid  Filled 
Projectiles,"  J.  of  Spacecraft,  Vol.  15,  No.  6,  348-354,  December  1978. 

4.  Kitchens,  C.W.,  Jr.  and  Gerber,  N.,  Prediction  of  Spin  Decay  of  Liquid  Filled 
Projectiles,  BRL  Report  1996,  Aberdeen  Proving  Ground,  Md.,  July  1977. 

5.  Kitchens,  C.W.,  Jr.,  Gerber,  N.,  and  Sedney,  R.,  Oscillations  of  a  Liquid  in  a 
Rotating  Cylinder:  Part  I.  Solid-Body  Rotation.  ARBRL-TR-02081,  Aberdeen 
Proving  Ground,  Md.,  June  1978. 

6.  Kitchens,  C.W.,  Jr.,  Gerber,  N.,  and  Sedney,  R.,  Oscillations  of  a  Liquid  in  a 
Rotating  Cylinder:  Part  II.  Spin-Up,  BRL  Report  (in  preparation). 

7.  McLeod,  A.R.,  "The  Unsteady  Motion  Produced  in  a  Uniformly  Rotating  Cylinder 
of  Water  by  a  Sudden  Change  in  the  Angular  Velocity  of  the  Boundary", 
Philosophical  Magazine  and  Journal  of  Science,  Vol.  44,  No.  259,  1-14,  1922. 
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that  viscous  diffusion  was  the  dominant  spin-up  mechanism,  and  so  neglected  the 
nonlinear  convective  terms  in  the  equations.  Winch2  later  extended  McLeod's  arguments 
to  the  case  where  a  free  surface  was  present.  In  both  cases  the  theory  significantly 
overestimated  the  fluid  spin-up  time.  Greenspan,  et  al.6>6  considered  the  situation  in 
which  Rs  »  aRx/-J  is  large,  which  is  typical  for  the  cases  of  interest  here.  Convective 
terms  were  retained,  but  in  the  limit  where  only  small  changes  in  rotational  speed 
were  allowed,  i.e.,  small  Rossby  number,  £  -  aa/a.  The  latter  served  as  the  expan¬ 
sion  parameter  for  a  small  perturbation  analysis  of  the  problem.  It  was  found  that 
the  Ekman  layers  on  the  endwalls  drive  a  strong  secondary  flow  which  accelerates  the 
spin-up  process.  The  characteristic  spin-up  time  was  shown  to  scale  as  ts~  £  fie  ftA 


Wedemeyer10  was  the  first  to  carry  out  a  nonlinear  analysis  of  this  problem 
which  did  not  hinge  on  the  Rossby  number  being  small.  He  argued  that  the  direct 
influence  of  viscosity  was  likely  to  be  felt  only  in  the  Ekman  layers.  The  bulk  of  the 
flow  in  the  interior,  or  core  flow  as  he  called  it,  was  thus  assumed  to  be  inviscid. 
Further,  the  radial  and  axial  flow  velocities,  u  and  w,  would  be  expected  to  be  much 
smaller  than  the  azimuthal  component,  v,  in  this  region.  He  also  assumed,  based  on 
Benton's11  work,  that  the  Ekman  layer  structure  would  accommodate  itself  in  a  quasi¬ 
steady  fashion  to  changes  in  the  core  flow.  With  these  approximations,  the  Navier- 
Stokes  equations  reduce  to  a  single  nonlinear  partial  differential  equation  for  v  as  a 
function  of  r  and  t,  which  has  a  simple  closed  form  solution  for  the  case  of  impulsive 
spin-up  from  rest. 


Wedemeyer's  phenomenological  analysis  formed  the  groundwork  for  several 
amplifications  and  extensions.12-16  Briefly,  the  qualitative  picture  that  emerges 


8.  Greenspan,  H.P.  and  Howard,  L.N.,  "On  a  Time-Dependent  Motion  of  a  Rotating 
Fluid,"  J.  of  Fluid  Mechanics,  Vol.  17  Pt.  3,  385-404,  1963. 

9.  Greenspan,  H.P.  and  Weinbaum,  S.,  "On  Non-Linear  Spin-Up  of  a  Rotating  Fluid," 
J.  of  Mathematics  and  Physics,  Vol.  44,  66-85,  1965. 

10.  Wedemeyer,  E.H.,  "The  Unsteady  Flow  Within  a  Spinning  Cylinder,"  J.  of  Fluid 
Mechanics,  Vol.  20  Pt.  3,  383-399,  1964. 

11.  Benton,  E.R.,  "On  the  Flow  Due  to  a  Rotating  Disk,"  J.  of  Fluid  Mechanics, 
Vol.  24  Pt.  4,  781-800,  1966. 

12.  Benton,  E.R.,  "Vorticity  Dynamics  in  Spin-Up  from  Rest”,  Physics  of  Fluids, 

Vol.  22  No.  7,  1250-1251,  1979. 

13.  Venezian,  G.,  "Non-Linear  Spin-Up",  Topics  in  Ocean  Engineering,  Vol.  2,  87-96, 
Gulf  Publishing  Co.,  Houston,  1970. 

14.  Weidman,  P.D.,  "On  the  Spin-Up  and  Spin-Down  of  a  Rotating  Fluid.  Part  I: 
Extending  the  Wedemeyer  Model.  Part  II:  Measurements  and  Stability",  J.  of 
Fluid  Mechanics,  Vol.  77  Pt.  4,  685-735,  1976. 
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is  as  follows.  Due  to  the  no-slip  boundary  condition,  the  fluid  motion  in  the  Ekman 
layer  closely  follows  that  of  the  wall  beneath  it.  Thus,  when  the  cylinder  is  impulsively 
started,  this  fluid  is  flung  out  toward  the  sidewall  by  centrifugal  forces;  the  radial 
pressure  gradient  needed  to  prevent  such  a  flux  has  not  yet  had  time  to  develop,  since  it 
is  determined  by  the  still  more  or  less  quiescent  core  flow.  Since  fluid  volume  must 
be  conserved,  this  outward  flux  in  the  Ekman  layer  must  be  balanced  by  an  influx  of 
fluid  from  the  core  flow.  The  Ekman  layer  thus  provides  a  pumping  action  on  the 
interior  flow.  Near  the  cylindrical  axis  fluid  is  sucked  into  the  viscous  layer  where  it 
acquires  additional  angular  momentum,  is  flung  outward,  and  then  turned  back  into  the 
core  flow  to  complete  the  circuit. 

Hence  a  strong  secondary  flow  is  set  up.  The  core  flow  is  characterized  by  a 
(suitably  normalized)  azimuthal  velocity  profile,  v,  of  order  one,  and  a  radial  velocity, 
u,  of  order  (R/H)Re-1/2,  which  is  directed  in  toward  the  axis.  The  strong  rotational 
effects  in  this  flow  allow  application  of  the  Taylor-Proudman  theorem,17  which  indicates 
that  u  and  v  are  both  independent  of  z  to  this  approximation.  On  the  other  hand,  the 
axial  velocity,  w,  varies  linearly  from  zero  at  the  mid-plane  to  a  maximum  value  at 
the  endwall,  where  it  is  accommodated  to  the  wall  through  the  Ekman  layer  structure. 
It  is  found  that  w  also  changes  sign  between  the  axis  and  the  sidewall  to  reflect  the 
sink/source  nature  of  the  matching  with  the  Ekman  layer. 

Thus,  while  viscosity  is  responsible  for  transmitting  angular  momentum  from  the 
cylinder  to  the  fluid  in  the  Ekman  layer,  it  is  the  inviscid  core  that  is  responsible  for 
redistributing  it  ova1  the  rest  of  the  flow.  It  should  be  noted  that  the  characteristic 
viscous  diffusion  time,  tj),  scales  as  0(kV*>)  =  Re  a ',  and  hence  IV 1 «  ts  «  t  .  The  fluid 
is  seen  to  achieve  solid-body  rotation  on  a  much  shorter  time  scale  than  tD;  hence  it 
really  is  the  secondary  flow  and  not  viscous  diffusion  alone  that  is  responsible  for  spin- 
up,  as  evidenced  by  the  poor  correlation  of  the  analyses  of  McLeod7  and  Winch2  with 
experiment. 


15.  Watkins,  W.B.  and  Hussey,  R.G.,  "Spin-Up  from  Rest:  Limitations  of  the 
Wedemeyer  Model",  Physics  of  Fluids,  Vol.  16  No.  9,  1530-1531,  1973. 

16.  Watkins,  W.B.  and  Hussey,  R.G.,  "Spin-Up  from  Rest  in  a  Cylinder",  Physics  of 
Fluids,  Vol.  20  No.  10,  1596-1604,  1977. 

17.  Greenspan,  H.P.,  The  Theory  of  Rotating  Fluids,  Cambridge  University  Press, 
London,  1968. 
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Recent  advances  in  computer  hardware  and  algorithms  have  made  it  feasible  to 
attack  this  problem  from  first  principles.  Thus  Warn-Varnas,  et  al.,18  and  Kitchens18>20 
have  each  solved  the  full  unsteady,  axisymmetric  Navier  Stokes  equations  for  the  spin- 
up  of  a  viscous  incompressible  fluid  in  a  completely  filled  cylinder,  using  finite-difference 
techniques.  Their  calculations  agree  with  one  another,  and  with  laboratory  experiments. 
They  further  reinforce  the  qualitative  picture  described  above  based  on  the  partly 
phenomenological  Wedemeyer  model. 

The  investigations  discussed  above,  with  the  exception  of  that  by  Winch, 2  were 
all  concerned  with  the  case  where  the  liquid  completely  fills  the  cavity.  In  practice, 
however,  the  shells  are  often  only  partially  filled.  Moreover,  field  detail  show  that 
the  probability  of  a  projectile  experiencing  an  erratic  flight  is  a  strong  function  of  its 
fill  ratio,  i.e.,  the  percentage  of  the  cavity  volume  occupied  by  fluid.  Thus  there  is 
an  important  need  to  understand  the  influence  of  the  interior  free  surface  on  the  spin- 
up  process. 

Gerber22’2^  has.  shown  that  the  final  form  of  this  surface,  after  solid-body 
rotation  is  achieved,  is  a  parabola  whose  shape  is  determined  by  the  fill  ratio  and  the 
Froude  number,  F  =  (aR)1/$  H  .  Winch's2  attempt  to  predict  the  evolution  of  this 
shape  with  time  failed  because  he  neglected  the  dominant  secondary  flow.  Goller  and 
Ranov24  have  extended  Wedemeyer's  model  to  include  the  free  surface  motion,  and 
were  able  to  show  that  it  acts  to  retard  the  spin-up  process.  They  presented  numerical 
results  for  both  the  free  surface  shape  and  the  azimuthal  velocity  profile  as  functions 


18.  Warn-Varnas,  A.,  Fowlis,  W.W.,  Piacsek,  S.,  and  Lee,  S.M.,  "Numerical  Solutions 

and  Laser-Doppler  Measurements  of  Spin-Up",  J.  of  Fluid  Mechanics,  Vol.  85 
Pt.  4,  609-639,  1978. 

19.  Kitchens,  C.W.,  Jr.,  Navier-Stokes  Solutions  for  Spin-Up  from  Rest  in  a 

Cylindrical  Container,  ARBRL-TR-02193,  Aberdeen  Proving  Ground,  Md., 
September  1979. 

20.  Kitchens,  C.W.,  Jr.,  Navier-Stokes  Solutions  for  Spin-Up  by  a  Predictor  Corrector 
Multiple-Iteration  Method,  AIAA  Paper  No.  79-1454,  presented  at  AIAA 
Computational  Fluid  Dynamics  Conference,  July  23-25,  1979,  Williamsburg,  Va. 

21.  Mark,  A.,  Measurements  of  Angular  Momentum  Transfer  in  Liquid-Filled 

Projectiles,  ARBRL-TR-2029,  Aberdeen  Proving  Ground,  Md.,  November  1977. 

22.  Gerber,  N.,  "Properties  of  Rigidly  Rotating  Liquids  in  Closed  Partially  Filled 
Cylinders",  ASME  Transactions,  J.  of  Applied  Mechanics,  Vol.  97,  734-735,  1975. 

23.  Gerber,  N.,  Rigidly  Rotating  Liquids  in  Closed  Partially  Filled  Cylindrical 

Cavities,  BRL-MR-2462,  Aberdeen  Proving  Ground,  Md.,  March  1975.  ~ 

24.  Goller,  H.  and  Ranov,  T.,  "Unsteady  Rotating  Flow  in  a  Cylinder  with  a  Free 

Surface",  ASME  Transactions,  J.  of  Basic  Engineering,  Vol.  90D,  No.  4, 
445-454,  December  1968. 
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of  time.  The  former  were  compared  with  data  from  a  series  of  laboratory  experiments 
and  th.e  agreement  was  quite  good.  The  problem  with  applying  their  analysis  to  actual 
configurations  is  that  Goller  and  Ranov  assumed  the  free  surface  did  not  intersect 
either  of  the  endwalls.  This  is  a  rather  restrictive  assumption,  as  shown  by  the  analysis 
of  Gerber.22,23  The  purpose  of  the  investigation  reported  here  is  to  extend  the 
Wedemeyer-Goller-Ranov  model  to  situations  where  the  free  surface  intersects  one  or 
both  of  the  endwalls. 

The  analysis  is  described  in  Section  2,  and  numerical  results  are  presented  and 
discussed  in  Section  3.  The  conclusions  drawn  are  summarized  in  Section  4.  To  help 
preserve  continuity  in  the  main  text,  the  equations  describing  the  free  surface  motion 
are  developed  in  Appendix  A  and  called  out  as  needed.  Appendix  B  provides  some 
details  on  the  finite  difference  algorithm  used  in  solving  the  equations. 


Section  2 
ANALYSIS 


We  wish  to  predict  the  spin-up  of  fluid  in  a  cylinder  of  radius  R  and  height  H 
(see  Figure  1)  which  at  t  =  0  impulsively  begins  to  spin  with  constant  angular  velocity, 
Jl  ,  about  its  axis.  The  fluid  is  incompressible  and  characterized  by  its  density,  p  , 
and  kinematic  viscosity,  V  .  Its  initial  level  when  at  rest  is  denoted  by  L.  The 
velocity  components  in  the  coordinates  are  (u,  v,  w)  respectively. 


The  fluid  motion  is  governed  by  the  axisym  metric  Navier-Stokes  equations. 
However,  following  Wedemeyer,!®  the  motion  in  the  core  flow  is  assumed  to  be 
predominantly  rotational  in  nature,  so  that  terms  involving  u  and  w  are  neglected.  This 
leads  to  a  Taylor-Proudman  type  flow,  i.e.,  to  this  degree  of  approximation  u,  v  and  p 
are  functions  of  only  r  and  t.  After  converting  to  dimensionless  variables  (see 
Nomenclature  list),  the  equation  for  v  (r,  t)  reduces  to, 


3_V 

3t 


+  UL 


/ii.  *  £.) 

\  dr  r  ) 


subject  to  the  initial  and  boundary  conditions, 

O  *  r  4  i 
t  >  0 
t  >  o 


V  (rt  O')  *0 

V  (_o ,  t )  =  0 

ir  (i ,  t  )  *  i 


(l) 


(2a) 

(2b) 

(2c) 


An  additional  expression  relating  u  to  v  is  required  to  close  the  system.  Assuming 
for  the  moment  that  this  is  known,  we  proceed  to  form  a  finite  difference  analog  to 
Eq.  (1),  using  a  uniform  grid  in  space  and  time.  Thus,  we  define  rj  =  (j-1)  A  r, 
j  =  1,2  .  .  .  NJ  where  A  r  =  1/(NJ-1).  Then  denoting  v  (rj,  t  =  i  At)  by  v'.  >  we  can 
express  the  solution  at  the  next  time  step  by 


V 


t 


V 1  +  ~  £G-4(u,i/)  «f  Ql  +  \u.,  ir  )  J  +  0(At)3 


(3a) 
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STAGE  1 


STAGE  3 


Figure  1  CYLINDRICAL  GEOMETRY  AND  POSSIBLE  FLUID 
CONFIGURATIONS  FOR  L/H  >  % 


where 


Centered  second-order  accurate  differences  are  used  to  approximate  the  spatial 
derivatives  in  Eq.  (3b).  When  assembled  into  Eq.  (3a),  the  right-hand  side  will  contain 
v^1 ,  v*1  ,  and  .  Thus  the  system  is  implicit  in  that  all  velocity  values  at  the 
new  time  step  must  be  solved  for  simultaneously.  It  is  also  tridiagonal,  and  so  lends 
itself  to  solution  by  the  efficient  Thomas  algorithm;25  the  details  are  given  in  Appendix 
B.  An  implicit  scheme  offers  the  advantage  of  using  time  steps  which  are  ordinarily 
larger  than  would  be  allowed  by  an  explicit  scheme,  such  as  that  used  by  Goller  and 
Ranov.24 


The  presence  of  u  in  Eq.  (3  b)  represents  a  nonlinearity,  which  is  handled  via  an 
iterative  relaxation  procedure.  The  distributions  v1?1  and  vf  are  presumed  known,  and 
by  linear  extrapolation  a  first  guess  is  made  for  v*j  *  .  From  this,  an  estimate  of  u1?1 
is  made,  and  then  Eq.  (3)  is  solved.  If  at  each  point  this  solution  agrees  with  that 
which  was  estimated  to  within  some  prescribed  relative  error,  it  is  accepted  and  all 
variables  are  updated  in  anticipation  of  the  next  time  step.  If  they  do  not  agree,  the 
solution  is  relaxed  as  follows 

(r*  ■*  1  )  S  ,  (rrO 

ir.  =  air.  +  ( 1  -  OC  )  V  ■  (5) 

4  i  1 

where  v^  represents  the  mth  iterate  to  v'J1,  and  v?  represents  the  solution  vector  to 
Eq.  (3)  based  on  this  iterate.  This  is  repeated  until  the  desired  convergence  is  achieved. 
It  was  found  necessary  to  vary  the  relaxation  parameter,  «_  ,  during  the  course  of 
the  calculation  to  insure  convergence.  The  values  used  are  given  in  Appendix  B. 


25. 


Richtmeyer,  R.D.  and  Morton,  K.W.,  Difference  Methods  for  Initial  Value 
Problems,  2nd  Ed.,  Chapter  I,  Interscience  Publishers,  New  Vork,  1967. 
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The  above  discussion  presumes  the  existence  of  a  relation  determining  u(r,  t)  as 
a  function  of  v(r,  t).  For  the  fully-filled  cylinder,  Wedemeyer  argued  that  to  preserve 
continuity  the  radially  outward  flux  in  the  endwall  Ekman  layers  must  be  balanced  at 
each  radius  by  an  equal  but  opposite  inward  flux  in  the  core  flow.  Further,  Benton11 
had  shown  that  the  Ekman  layer  flow  established  itself  on  a  time  scale  of  the  order  of 
a  few  revolutions  of  the  cylinder,  which  is  orders  of  magnitude  faster  than  the  interior 
fluid  spin-up  time.  Hence  the  efflux  in  the  layers  was  assumed  to  accommodate  itself  in 
a  quasi-steady  manner  to  changes  in  the  core  flow.  Due  to  the  columnar  nature  of 
the  motion  the  radial  influx  in  the  core  is  spread  uniformly  with  depth.  This  led 
Wedemeyer  to  the  following,  as  he  called  it,  compatibility  condition  relating  u(r,  t)  to 
v(r,  t): 

ul  =  -O.88<0  [-jj-J  Re  (r-lr)  (6) 


For  a  partially-filled  cylinder,  there  is  an  additional  contribution  to  u  resulting 
from  the  free,  surface  motion.  The  drop  in  fluid  level  near  the  axis  and  the  accompanying 
rise  at  the  sidewall  imply  a  radial  outflux  in  opposition  to  the  influx  generated  by  the 
Ekman  layers.  Further,  these  layers  will  no  longer  necessarily  cover  the  entire  endwall. 
These  features  of  the  flow  serve  to  retard  the  spin-up  process.  For  convenience  we 
separate  the  contributions  to  u  from  the  top  and  bottom  Ekman  layers  and  the  free 
surface  motion  as  follows, 


U~ 


+  ^BEL  +  U  FS 


(7) 


The  form  each  of  these  will  take  at  any  instant  depends  on  whether  the  free 
surface  intersects  neither,  the  top,  or  both  of  the  endwalls.  Following  Gerber, 22, 23 
we  shall  refer  to  these  successive  configurations  as  Stages  1,  2b,  and  3  as  indicated 
in  Figure  1.*  The  remaining  discussion  naturally  divides  itself  along  the  same  lines. 


Gerber  also  treated  a  fourth  configuration,  Stage  2a,  in  which  the  free  surface 
intersects  the  bottom,  but  not  the  top,  endwall.  This  occurs  only  for  fill  ratios, 
L/H,  less  than  one-half.  Such  cases  are  not  considered  here,  although  the 
extension  of  the  present  methods  to  flows  in  Stage  2a  should  be  straightforward. 
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Stage  1 


This  is  the  ease  treated  previously  by  Goller  and  Ranov,24  in  which  neither 
endwall  is  intersected.  Before  u  can  be  found,  one  must  first  determine  the  free 
surface  contour,  Zpg(r,t),  given  v(r,t).  This  is  done  using  Eqs.  (A-3)  and  (A-5)  from 
Appendix  A, 


Z  At)  = 

L_ 

+  F  f  ir*  r  dr 

(8a) 

i 

H 

< 

3 

Z  -  F  f  ~  dr' 
i  J  r' 

(8b) 

r 


Simpson’s  rule  is  used  to  evaluate  the  quadatures  on  the  radial  grid. 


Clearly  for  this  stage,  u^el  =  0.  The  bottom  Ekman  layer,  however,  is  assumed 
to  pump  fluid  at  the  same  rate  as  for  the  fully  filled  case.  But  whereas  before  the 
compensating  flux  in  the  core  flow  was  spread  uniformly  over  half  the  cylinder  height 
(by  symmetry),  it  is  now  spread  over  the  fluid  depth,  Zpg: 


27 Tr  2 


fs  ubcl 


Z 


6 

u  cLz 


(9) 


For  the  Reynolds  numbers  of  interest  here,  the  Ekman  layer  flow  is  expected  to 
be  laminar;  for  such  cases  the  flux  integral  on  the  right  has  previously  been  represented 
approximately  as  one  or  the  other  of  the  following 


1} 


c Li 


?  c %>a)v* 


-  i 


0-4-+  3  O  ~  ui) 


1  '  n 

—  Z 

*  *  *o 


(10a) 

(10b) 


The  upper,  linear  expression  in  Eq.  (10a)  was  used  by  Wedemeyer.1®  The  series  given 
in  Eq.  (10b)  is  a  polynominal  fit  by  Goller  and  Ranov24  to  the  more  refined  numerical 
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data  of  Rogers  and  Lance. 26  An  analogous  expression  for  this  flux  has  also  been 
derived  for  situations  where  the  Ekman  layer  flow  is  expected  to  be  turbulent.10 
Substituting  Eq.  (10)  into  (9)  gives,  in  dimensionless  form, 


“«*.  -  -a-  ir  /“->■  (11) 

n  i- 

As  for  the  free  surface  motion,  the  flux  across  any  radial  surface,  r,  should 
equal  the  rate  of  change  of  fluid  volume  between  the  axis  and  7.  Again,  this  is  spread 
uniformly  over  the  height  of  the  fluid  column: 


Z  It  r 


r 


/ 


or  in  dimensionless  form. 


u 


FS 


(12) 


where  backward  differencing  is  used  for  dZ^/dt, 


o  7 

Izfs.) 

A 


7  _  2  c 

A  t 


+  O(At) 


(13) 


and  Simpson's  rule  is  again  used  for  the  quadrature.  Use  of  Eqs.  (11)  -  (13)  in  (7) 
completes  the  evaluation  of  u  in  terms  of  v  needed  in  Eq.  (3b).  The  above  analysis 
of  Stage  1  conditions  is  identical  to  that  used  by  Goller  and  Ranov,  except  for  the  use 
of  an  implicit,  as  opposed  to  their  explicit,  time-marching  scheme. 


It  is  shown  in  Appendix  A  that  for  F  >  4(1  -  L/H)  the  flow  will  at  some  time 
intersect  the  top  endwall.  This  will  be  evidenced  by  Eq.  (8a)  predicting  Zf  >  l.  The  flow 
is  now  in  Stage  2b,  which  is  discussed  next. 


26.  Rogers,  M.H.,  and  Lance,  G.N.,  "The  Rotationally  Symmetric  Flow  of  a  Viscous 
Fluid  in  the  Presence  of  an  Infinite  Rotating  Disk",  J.  of  Fluid  Mechanics, 
Vol.  7  Pt.  4,  617-631,  1960. 


Stage  2b 

Suppose  that  the  point  at  which  the  integral  in  Eq.  (8a)  first  exceeds  (1  -  L/H) 
lies  in  the  interval  between  grid  points  rj  =  rL  and  rj+i  =  rR,  i.e.,  pL  <  rR  *  rR. 
Using  linear  interpolation,  the  integrand  is  expressible  as 

^  ( r )  h  yz  r  -  b  +  2  a.  C  r  -  rL)  &  r  ±  rR  (14) 

where 

*  ■  /(»*r) 

If  Eq.  (14)  is  substituted  into  Eq.  (8a),  and  Zj  is  set  to  unity,  the  following  quadratic 
is  obtained  for  rH: 


a-F  (  rH  -rL)*  +  bF  ( rH-ru )  +  F1L  -  (i  -  L/H  )  =  o 

(  rL 

where  Zu  =  J  $(r)  dr  .  The  appropriate  root  (rjj  >  r^)  is 

,  _  „  .  -bFt  {(bn*  -<K*.F)[FI.-(I-  L/H)]y/t 

rM  ru  +  - - - - ■ —  'I5) 

2.  a.  F 

This  represents  the  radius  of  intersection  with  the  top  endwall.  For  r  >  rjj,  Zpg  =  1 
(see  Figure  1)  while  for  r  <  in  Eq.  (A-ll)  is  used, 

2„(r,t)  .  i-F  f"  S.'  c ir'  (16) 

r  r 

integrating  in  towards  the  axis,  r  =  0. 


Applying  the  same  arguments  as  in  Stage  1, 

,r  Si 


IX 


F  S 


(ri 


we  find  for  upg 

ps  )  OS  r  4  rH 


(17a) 


u*w%  ■  o  rH  4  r  4  1  (17b) 

One  can  show  from  continuity  arguments  that  the  integral  in  Eq.  (17a)  vanishes  at 
r  =  rjj»  and  hence  upg  remains  continuous  there. 
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Next  consider  the  Ekman  layer  pumping.  The  contribution  from  the  bottom  wall 
is  the  same  as  in  Stage  1, 

u  =  -Rc~'‘  IT  T  <18) 

PS 

But  the  top  wall  is  only  partially  wetted.  The  author  is  unaware  of  any  published 
investigations  into  the  nature  of  the  boundary  layer  that  forms  in  such  a  case.  At 
the  very  least,  one  expects  that  the  outward  flux  in  the  layer  must  vanish  at  ry,  since 
there  is  no  longer  any  fluid  at  r  <  r^  to  sustain  it.  So  to  the  mass  flux  in  Eq.  (10), 
and  hence  is  added  a  pre-multiplicative  factor,  KT  T  ,  to  reflect  this: 

u._._  ~  o  o  4  r  4  ru  (19a) 

T  cL  n 

urri  =  -  Re  >z  -&■  r  Kt  (r-rH)*T  f  (oj)  rH  &  r  4  i  (19b) 

where  KT,  7fT  >  0. 

Such  a  heuristic  extension  of  the  arguments  put  forth  by  Wedemeyer  and  Goller 
and  Ranov  satisfies  the  following  constraints: 

•  u  remains  continuous  across  r  =  rjjj 

•  the  contribution  from  u^el  vanishes  at  both  r  =  rjj  and  1,  as  it  should; 

•  the  influence  of  the  top  endwall  Ekman  layer  will  increase  as  rj]  decreases, 
which  is  intuitively  satisfying; 

•  the  columnar  nature  of  the  flow  is  maintained. 

Certainly  other  analytical  forms  can  be  found  which  would  satisfy  these  constraints. 
But  in  the  absence  of  any  experimental  data  to  serve  as  a  guide,  it  was  decided  to 
employ  the  simplest  form  possible.  For  the  calculations  reported  here,  KT  and 
were  both  set  to  unity. 


Eq.  (A-28)  indicates  that  for  F  >  (1  -  L/H)"*  eventually  the  fluid  displacement 
will  expose  a  portion  of  the  bottom  endwall.  This  will  be  indicated  during  the  evaluation 
of  Eq.  (16)  if  at  some  r,  Zpg  becomes  negative.  The  flow  is  now  in  its  final  configuration, 
Stage  3. 

Stage  3 

Once  the  bottom  endwall  is  intersected,  one  cannot  expect  the  boundary  condition 
in  Eq.  (4a)  to  apply,  as  there  is  no  longer  any  fluid  at  the  axis.  Instead,  since  the 
fluid  depth  goes  to  zero  at  the  intersection  point,  r0,  it  seems  appropriate  to  impose  a 
no-slip  boundary  condition  there, 


v  (r0)  *  r0 


(20) 


and  to  include  only  those  grid  points  rj  >  r0  in  the  solution  of  Eq.  (3).  These  changes 
can  be  handled  with  a  straightforward  modi?4  cation  to  the  solution  algorithm  used  in 
Stages  1  and  2b  (see  Appendix  B). 

A  description  of  the  free  surface  contour  as  z  =  Zpg(r,  t)  can  be  developed  for 
this  stage  just  as  was  done  in  the  ear'ier  Stages.  But  such  a  development  suffers  from 
two  significant  disadvantages.  Whereas  during  Stage  2b  we  had  only  one  time-varying 
limit  of  integration,  rH,  during  Stage  3  we  have  the  added  complication  of  r0  also 
changing  with  time.  Moreover,  it  is  not  possible  to  represent  either  of  these  limits 
explicitly  in  terms  of  the  given  velocity  profile.  Secondly,  the  free  surface  in  this  stage 
is  typically  quite  steep,  so  that  the  accuracy  of  the  numerical  quadratures  on  the  radial 
grid  comes  into  question. 


Both  obstacles  may  be  overcome  by  turning  our  frame  of  reference  sideways. 
That  is,  we  can  just  as  well  choose  z  as  our  independent  variable  and  let  r  =  Rp§(z,  t) 
represent  the  surface  contour  (see  Figure  1).  In  this  frame  the  fluid  is  always  constrained 
between  the  constant  limits  z  =  0  and  1.  More  important,  the  surface  has  a  shallow 
dope,  allowing  the  application  of  standard  quadrature  techniques.  From  Eqs.  (A-22) 
and  (A-24),  the  free  surface  will  now  evolve  according  to 
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(21a) 
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*  dz' 

Zr  2  (Z  t  ) 


(21b) 


where  v(z,  t)  is  understood  to  mean  v(Rpg(z,  t),  t).  This  defines  the  contour  implicitly. 

o 

For  a  given  v(r,  t),  Eq.  (21)  may  be  viewed  as  an  integral  equation  for  Rpg.  It  is  also 
weakly  nonlinear  owing  to  the  implicit  dependence  of  v(z,  t)  on  the  solution.  As  a 
result,  during  Stage  3,  in  addition  to  the  "outer"  iteration  on  v  used  to  solve  Eq.  (3), 
we  will  now  have  an  "inner"  iteration  process  to  determine  the  free  surface  contour 
consistent  with  v. 


To  solve  Eq.  (21),  a  uniform  grid  is  laid  out  in  the  z  direction,  zk  =  (k  -  1)  az, 
k  =  1,2  .  .  .  NK  where  Az  =  1/(NK  -  1).  On  the  initial  transition  to  Stage  3,  a  first 
guess  at  Rpg  is  made  by  fitting  a  cubic  spline27  through  the  ordered  pairs  (Zpgj,  rj) 
obtained  from  Eq.  (16).  This  spline  is  then  evaluated  at  all  the  zk,  which  specifies 
Rpsk=  Rps(zk,  t).  On  succeeding  steps  in  Stage  3,  the  first  guess  to  Rpg  is  taken  to 
be  that  from  the  previous  time  step.  Next,  a  cubic  spline  is  passed  through  the  point 
(r0>  f0>  (see  Eq*  (20))  and  the  ordered  pairs  (rj,  Vj);  again  only  those  rj  >  r0  are  included 
in  the  fit.  This  spline  is  then  evaluated  at  all  the  Rpsk>  which  defines  vk  =  v(zk,  t). 

Having  thus  specified  Rpg(z,  t)  and  a  self-consistent  v(z,  t)  in  discrete  form, 
the  quadratures  in  Eq.  (21)  are  performed  using  Simpson's  rule.  Convergence  is  achieved 
when  this  solution  agrees,  at  each  k,  with  the  previous  estimate  to  within  some  prescribed 
relative  error.  If  not,  the  solution  is  relaxed  according  to 

(n  +  t)  «  (n) 

3;/?«  (22) 

(n)  s 

where  Rpg  represents  the  nth  iterate,  and  Rpg  represents  the  solution  vector  to  Eq. 
(21)  based  on  this  iterate.  Note  that  to  remain  self-consistent  the  velocity  spline 
defining  ^  must  be  refit  and  reevaluated  at  each  iteration.  For  linear  equations,  it 
can  be  shown  (e.g.,  Reference  28)  that  ft  =  1  will  yield  converged  solutions.  In  our 
case,  the  nonlinearity  necessitated  underrelaxing  the  iterations;  a  value  of  =  0.5  was 
used  and  typically  no  more  than  4-5  iterations  were  required  for  convergence.  When 
a  solution  is  reached,  cubic  splines  are  again  employed  to  revert  from  Rpg(zk,  t)  back 
to  Zpg(rj,  t). 

27.  IMSL  Library  Reference  Manual,  8th  Ed.,  Vol.  1,  Chapter  I,  IMSL  Inc.,  Houston, 
Texas,  1980. 

28.  Hildebrand,  F.B.,  Methods  of  Applied  Mathematics,  2nd.  Ed.,  Chapter  3,  Prentice- 

Hall,  Inc.,  Englewood  Cliffs,  N.J.,  1965^ 
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Using  the  same  physical  arguments  as  in  the  earlier  stages,  but  with  r  replaced  by 
z  as  the  integration  variable,  leads  to  the  following  equation  for  upg: 


ups  '  / 


'n  3  R 


FS 


PS 


dt 


cLz!  / (riF s) 


where,  analogous  to  Eq.  (13), 

i  +  1 


3  R, 
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dt  Jk 
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q  *  1  n  <• 
P-FSf  ~  *P 


At 


C  ^  ^  *  rH 

(23a) 

rH  *  r  4  i 

(23b) 

+  O  (At) 

(24) 

It  is  clear  that  upg  vanishes  at  r  =  r0  (Zpg  =  0);  one  can  also  show  from  conservation 
of  fluid  volume  that  upg  =  0  at  r  =  rg  (Zpg  =  1),  as  it  should.  In  evaluating  Eq.  (23), 
first  Simpson's  rule  is  used  to  perform  the  integration  with  the  upper  limit  replaced 
by  zjj,  for  all  k  =  1,2  .  .  .  NK.  A  cubic  spline  is  then  fit  through  these  flux  integrals, 
effectively  defining  the  flux  as  a  continuous  function  of  z.  The  value  of  this  spline 
at  Zpgj  then  represents  the  numerator  in  Eq.  (23)  needed  to  obtain  upg^.  It  should 
be  emphasized  here  that  the  axial  grid  is  used  only  to  evaluate  Eqs.  (21)  and  (23) 
pertaining  to  the  free  surface  motion.  Eq.  (3)  for  the  azimuthal  velocity  is  still  solved 
on  the  radial  grid. 

The  radial  velocity  contribution  from  the  top  Ekman  layer  is  still  given  by  Eq. 
(19).  The  contribution  from  the  bottom  layer  is  now  also  presumed  reduced  by  a  pre- 
multiplicative  factor  to  reflect  the  fact  that  it  no  longer  covers  the  entire  wall: 

U-Bel  ■  °  r  i  re  (25a) 

U-Bml  m  '  Re'  ^  7 T  (  1  *  T  )*  f >  ro^r^1  (25b) 

where  Kg,  /g  >  0;  both  were  set  to  unity  for  the  calculations  presented  here.  This 
expression  satisfies  our  intuitive  expectations  that  ugpL  should  vanish  at  r0,  and  that 
its  contribution  at  any  r  should  grow  progressively  smaller  as  more  of  the  wall  is 
exposed. 
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In  closing  this  section,  we  wish  to  recognize  that  the  model  described  above, 
heuristic  as  it  is,  forms  a  rather  simplified  picture  of  a  complex  flowfield.  But  the 
equations  are  felt  to  adequately  describe  the  gross  features  of  the  phenomena,  and  as 
such  should  provide  a  useful  first  step  towards  a  deeper  understanding.  We  are  buoyed 
in  this  belief  by  the  successful  correlations  with  experiment  of  the  earlier  simplified 
models  of  Wedemeyer10  and  Goller  and  Ranov.24 


Section  3 

NUMERICAL  RESULTS  AND  DISCUSSION 


To  uniquely  define  the  problem,  one  must  specify  values  for  Re,  F,  H/R  and 
L/H.  Values  of  these  input  parameters  for  the  cases  to  be  discussed  here  are  summarized 
in  Table  1,  along  with  the  temporal  and  radial  grid  sizes  used.  Representative 
calculations  were  made  for  smaller  values  of  A  r  and  At  to  insure  that  convergence 
had  been  achieved.  The  cases  include  flows  whose  final  configurations  lie  in  each  of 
the  three  permissible  stages.  For  those  cases  where  Stage  3  was  achieved  NK  =  51 
axial  grid  points  were  used  for  the  free  surface  integrations.  In  all  cases  Eq.  (10b) 
was  used  to  represent  the  nominal  Ekman  layer  flux.  A  relative  error  criterion  of 
10-3  was  used  for  both  the  velocity  iterations,  and  the  free  surface  iterations  during 
Stage  3. 


Table  1 

INPUT  PARAMETERS  FOR  CASES  1-6 


Case 

No. 

Re 

F 

H/R 

L/H 

Final 

Configuration 

At 

HI 

n 

2.293  x  105 

0.6746 

2.865 

0.8 

Stage  1 

9.251 

0.02 

H 

1.172  x  105 

4.141 

3.0 

0.8 

Stage  2b 

5.0 

0.02 

H 

1.172  x  105 

3.5 

3.0 

0.6 

Stage  3 

5.0 

0.02 

■ 

1.172  x  104 

3.5 

3.0 

0.6 

Stage  3 

2.0 

0.02 

5 

1.172  x  103 

3.5 

3.0 

0.6 

Stage  3 

0.5 

0.02 

6 

1.172  x  105 

6.0 

3.0 

0.8 

Stage  3 

5.0 

0.01 

A  problem  occurs  in  defining  when  the  fluid  is  fully  spun  up  which  is  similar  to 
that  which  arises  in  defining  the  "edge"  of  a  boundary  layer  -  the  condition  is  approached 
only  asymptotically.  As  shown  in  Appendix  A,  the  final  fluid  configuration,  and 
corresponding  angular  momentum,  depend  only  on  L/H  and  F.  Hence  the  latter  can 
be  computed  a  priori.  At  each  step  in  the  calculation,  the  current  angular  momentum 
is  evaluated  and  referenced  to  this  final  value.  The  instant  at  which  this  ratio  first 
equals  or  exceeds  0.99  will  be  used  here  as  the  definition  of  the  spin-up  time,  t§. 
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The  conditions  for  Case  1  were  chosen  to  match  those  given  in  Figures  10  and  14 
of  Reference  24,  which  include  experimental  data.  As  noted  earlier,  for  such  Stage  1 
flows  the  present  model  is  virtually  identical  to  Goller  and  Ranov's,  save  for  the  use 
here  of  an  implicit  algorithm.  The  odd  time  increment  in  Table  1  was  chosen  as  a 
convenient  sub-multiple  of  the  dimensional  time  increments  at  which  experimental  data 
were  taken.  The  predicted  azimuthal  velocity  profiles  and  free  surface  contours  are 
shown  in  Figures  2a  and  2b,  respectively,  at  selected  instants  of  time.  Note  in  Figure 
2b  that  the  axial  coordinate  has  been  renormalized  to  the  same  scale  as  the  radius  to 
better  reflect  the  true  picture.  The  dots  represent  data  read  from  Goller  and  Ranov's 
curves  for  t  =  100  and  400  sec.;  the  agreement  with  the  present  calculations  is  quite 
good.  No  experimental  data  were  taken  for  the  velocity  profiles.  The  theory  predicts 
that  solid-body  rotation  is  achieved  at  tg  =  570  sec.,  the  last  instant  shown. 

The  conditions  for  Case  2  were  chosen  so  that  the  final  fluid  configuration  would 
be  in  Stage  2b.  Velocity  and  free  surface  profiles  for  this  case  are  shown  in  Figure  3 
for  dimensionless  times  out  to  ts  =  4000  in  increments  of  400.  The  upper  endwall  is 
intersected  early  on  in  the  calculation,  at  t  =  200.  The  asymptotic  approach  to  solid- 
body  rotation  is  clearly  evident  in  the  bunching  up  of  the  curves  toward  the  end  of 
the  calculation.  Unfortunately,  no  quantitative  experimental  data  are  available  for 
flows  in  Stage  2b. 


For  the  linearized  spin-up  of  fluid  without  a  free  surface,  the  basic  time 
dependence  has  been  shown  to  be  exponential.**  The  question  was  raised  in  the  published 
discussion  following  Reference  24  of  whether  similar  behavior  would  be  found  for  cases 
such  as  those  studied  here,  where  nonlinear  convective  terms  are  retained  and  a  free 
surface  is  present.  If  such  is  the  case,  one  might  expect  the  fluid  angular  momentum 
to  evolve  according  to 


Lm(t)  -at/*.'* 

— = —  «  i  -  e 

Z-g  (e°) 


(26) 


where  "a"  is  a  dimensionless  constant.  The  denominator  of  the  exponential  is  suggested 
by  the  finding  that  for  linearized  flows,  tg  scales  with  Re1/2  (Reference  8).  If  the 
present  calculations  do  indeed  reflect  such  behavior,  then  a  semi-logarithmic  plot  of 
the  momentum  deficit,  l-Lz(t)Lz(<*>),  vs.  Re-*/2  t  should  yield  a  straight  line.  Such  a 
plot  is  presented  in  Figure  4  for  Cases  1  and  2.  The  symbols  represent  numerical  data 
from  the  present  program;  a  straight  line  has  been  faired  through  the  data  for  each 
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Figure  2a  AZIMUTHAL  VELOCITY  PROFILES  FOR  CASE  1 
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case.  It  is  somewhat  remarkable  to  see  the  degree  to  which  the  data  do  in  fact  adhere 
to  a  pure  exponential  behavior.  Evidently,  however,  the  Re^/2  time  scaling  is  not 
perfect,  as  the  two  cases  fail  to  collapse  to  a  single  line;  also,  F  and  H/R  were 
different  in  the  two  calculations. 

Cases  3  through  5  are  for  conditions  such  that  the  final  fluid  configuration  is 
in  Stage  3.  They  differ  only  in  the  values  specified  for  Re,  which  range  over  two  orders 
of  magnitude.  Figure  5  shows  the  velocity  and  free  surface  profiles  for  Case  3;  results 
for  the  other  two  are  qualitatively  similar.  The  first  nine  curves  in  Figure  5a  represent 
v(r)  at  dimensionless  times  t  =  400  through  3600  in  increments  of  400.  The  last  curve  is 
at  tg  =  4590;  Cases  4  and  5  produced  values  of  tg  =  1200  and  231,  respectively.  The 
last  five  curves  end  abruptly  at  r  =  r0  because  there  is  no  fluid  to  the  left  of  this 
point,  and  so  v  is  undefined  there. 

The  free  surface  profiles  in  Figure  5b  are  given  for  the  same  times.  Again  the 
asymptotic  nature  of  the  approach  to  solid-body  rotation  is  evident  in  both  figures.  It 
is  also  clear  that  the  surface  contours  become  quite  steep,  and  thus  warrant  the  extra 
effort  spent  in  Stage  3  in  switching  from  a  radial  to  an  axial  grid  for  the  surface 
integrals.  Again,  no  experimental  data  for  flows  in  this  stage  are  available  for 
comparison. 

An  anomaly  was  noted  in  the  radial  velocity  distribution  during  Stage  3  of  Cases 
3-5.  Ordinarily  u(r),  as  predicted  by  Eqs.  (7),  (19),  (23)  and  (25),  is  negative,  i.e.,  there  is 
a  net  inward  flux  in  the  core  flow.  However,  for  some  inboard  radii  it  was  noted  that 
u  became  positive.  Such  a  sign  change  is  inconsistent  with  a  columnar  flow  model  in 
the  absence  of  sources  or  sinks.  However,  rather  than  representing  any  fundamental 
flaw  in  the  numerical  algorithm,  such  behavior  merely  reflects  our  inadequate  knowledge 
of  the  magnitude  of  the  Ekman  pumping  embodied  in  Eqs.  (19)  and  (25).  Indeed,  in 
the  absence  of  quantitative  data  for  such  flows,  it  is  difficult  to  assess  the  degree 
to  which  the  columnar  flow  approximation  remains  valid  in  the  vicinity  of  the  intersection 
point  r0.  To  keep  this  from  further  influencing  the  solution,  at  such  points  u  was  set 
to  zero.  Since  by  this  time  the  fluid  had  acquired  in  excess  of  90%  of  its  final 
momentum,  it  is  doubtful  that  this  would  have  a  serious  effect  on  the  results. 

The  momentum  deficit  vs.  time  for  Cases  3-5  is  shown  in  Figure  6.  Two  flagged 
symbols  are  shown  for  each  case;  these  represent  the  instants  just  after  the  transition 
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Figure  5a  AZIMUTHAL  VELOCITY  PROFILES  FOR  CASE  3 
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Figure  6  ANGULAR  MOMENTUM  DEFICIT  VS.  TIME  FOR  CASES  3  THRU  5 


between  stages.  Note  that  the  numerical  data  for  Cases  3  and  4  each  follow  a  line 
of  constant  slope  during  the  first  two  stages,  as  was  the  case  in  Figure  4.  After  the 
transition  to  Stage  3,  however,  the  results  follow  a  line  with  a  reduced  slope,  indicating 
a  slower  rate  of  spin-up.  The  fact  that  the  first  transition  apparently  has  little  effect 
on  the  rate  of  spin-up  can  be  explained  by  noting  that  the  radial  extent  of  the  fluid 
layer  at  the  top  wall  in  Figure  5b  is  relatively  thin.  Hence  its  contribution  to  the 
Ekman  pumping  remains  small  compared  with  that  fron.  the  bottom  layer,  which  still 
completely  covers  the  wall.  Upon  making  the  transition  to  Stage  3,  however,  the 
strength  of  the  bottom  Ekman  layer  pumping  is  progressively  reduced  as  the  wall  is 
exposed,  and  a  slower  spin-up  rate  results. 

The  results  for  Case  5,  the  lowest  Re,  would  seem  to  represent  an  exception  to 
the  rule.  All  the  numerical  data  for  this  case  appear  to  follow  a  single  straight  line, 
with  no  kink  at  the  second  transition.  The  following  qualitative  argument  is  offered 
to  explain  this  behavior.  The  Reynolds  number,  broadly  speaking,  represents  the  ratio 
of  convective  to  viscous  forces.  For  Cases  3  and  4  Re  was  large  enough  for  the 
convective,  i.e.,  secondary  flow,  effects  to  dominate.  Since  the  secondary  flow  is 
driven  by  the  endwall  Ekman  layers,  it  is  not  surprising  that  the  spin-up  rate  should 
depend  on  which  stage  the  flow  is  in.  For  Case  5,  however,  Re  takes  on  its  smallest 
value,  and  evidently  viscous  diffusion  from  the  sidewall  plays  a  greater  role.  Clearly, 
for  this  mechanism  the  degree  to  which  the  endwalls  are  exposed  is  not  as  critical, 
and  consequently  the  spin-up  rate  is  more  uniform.  Finally,  it  is  interesting  to  note 
in  Figure  6  that  the  angular  momentum  ratios  at  the  times  at  which  transitions  occur 
are  nearly  independent  of  Re. 

Case  6  is  similar  to  Case  3  except  that  the  values  of  both  L/H  and  F  have  been 
raised.  The  velocity  and  free  surface  profiles  for  this  case  are  shown  in  Figure  7. 
The  first  nine  curves  represent  t  =  400  through  3600  in  increments  of  400,  while  the  last 
is  at  t§  =  4060.  As  expected,  the  higher  Froude  number  in  this  case  results  in  a 
steeper  surface  contour  than  was  exhibited  for  Case  3.  The  momentum  deficit  vs.  time 
is  plotted  in  Figure  8,  where  the  data  for  Case  3  have  also  been  repeated  for  comparison. 
During  the  first  two  stages,  both  cases  collapse  to  nearly  the  same  straight  line.  Note, 
however,  that  in  Case  6  the  free  surface  intersects  the  top  wall  much  sooner,  and  the 
bottom  wall  much  later,  reflecting  the  higher  fill  ratio.  A  discontinuity  in  slope  is 
again  exhibited  at  the  second  transition. 
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Figure  7a  AZIMUTHAL  VELOCITY  PROFILES  FOR  CASE  6 
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Figure  7b  FREE  SURFACE  CONTOURS  FOR  CASE  6 


The  run  times  for  these  calculations  varied  from  0.73  min.  for  Case  1  to  3.3 
mins,  for  Case  6,  using  an  IBM  370/3031  system.  But  such  numbers  are  to  a  certain 
extent  arbitrary,  depending  as  they  do  on  what  temporal  and  spatial  grid  sizes  were 
used.  A  more  meaningful  measure  is  the  CPU  time  per  time  step  per  grid  point.  This 
was  approximately  1.5  x  10"3  sec.  for  the  first  two  cases,  and  about  double  that  value 
for  the  remainder.  The  latter  no  doubt  took  longer  because  of  the  additional  iterations 
needed  to  define  the  free  surface. 
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Section  4 

SUMMARY  AND  CONCLUSIONS 


A  theoretical  analysis  is  presented  for  the  axisym metric  spin-up  of  fluid  in  a 
partially-filled  cylindrical  cavity.  The  analysis  represents  an  extension  of  the  earlier 
treatments  by  Wedemeyer,1®  and  Goller  and  Ranov,24  to  those  cases  where  the  liquid 
free  surface  may  intersect  one  or  both  endwalls.  Earlier  estimates  of  the  Ekman 
layer  pumping  of  the  secondary  flow  are  modified  heuristically  for  situations  where  the 
layer  no  longer  covers  the  entire  wall.  Also,  due  to  the  very  steep  free  surface 
contour  in  the  latter  stages  of  spin-up,  it  was  found  advantageous  to  develop  the  free 
surface  equations  in  an  axial,  rather  than  radial,  coordinate  frame. 

A  computer  program  has  been  developed  which  solves  the  governing  nonlinear 
equations  using  a  straightforward  finite-difference  algorithm.  The  code  predicts  both 
the  azimuthal  velocity  distribution  with  radius,  and  the  free  surface  contour,  as  functions 
of  time.  Predictions  for  the  latter  agree  with  Gerber's22  steady-state  analysis  as  the 
limit  of  solid-body  rotation  is  approached.  The  present  results  also  exhibit  very  good 
agreement  with  the  time-resolved  experimental  data  for  the  surface  contour  taken  by 
Goller  and  Ranov.24  Their  data  were  taken  only  for  situations  where  the  free  surface 
did  not  touch  the  endwalls.  At  present,  there  are  no  quantitative  data  against  which 
the  present  analysis  can  be  compared  for  cases  where  one  or  both  endwalls  are 
intersected. 

Nevertheless,  the  following  qualitative  conclusions  have  been  drawn  from  the 
theory.  Plotting  the  fluid  angular  momentum  deficit  vs.  Re1/2  t  appears  to  correlate 
the  numerical  data  reasonably  well,  although  the  collapse  is  by  no  means  perfect.  Such 
plots  indicate  that  the  angular  momentum  transfer  follows  a  simple  exponential  behavior 
in  time.  For  Re  in  the  range  104  -  105,  the  growth  rate  appears  uniform  up  to  the 
point  where  the  bottom  endwall  is  intersected.  After  this,  exponential  behavior  is 
still  exhibited,  but  at  a  reduced  rate,  reflecting  the  diminishing  influence  of  the  bottom 
Ekman  layer  which  is  primarily  responsible  for  the  secondary  flow.  For  Re  on  the  order 
of  103,  a  uniform  growth  in  momentum  is  exhibited  over  the  whole  of  the  calculation, 
perhaps  suggesting  that  such  flows  are  dominated  by  viscous  diffusion  from  the  sidewall. 

To  validate  the  analysis  in  a  quantitative  sense,  it  is  recommended  that  its 
predictions  be  compared  with  suitably  designed  laboratory  experiments  and/or  more 
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refined  numerical  calculations,  e.g.,  solutions  of  the  full  axisymmetrie  Navier-Stokes 
equations.  Such  comparisons  could  also  serve  as  guidelines  for  adjusting  the  empirical 
parameters  in  Eqs.  (19)  and  (25)  used  to  define  the  strength  of  the  Ekman  pumping.  It 
is  felt  that  the  present  computer  code  could  then  serve  as  an  efficient  means  of 
predicting  fluid  spin-up  in  cases  of  practical  interest.  Such  predictions  form  an  integral 
part  of  any  attempt  to  predict  the  spin-decay  of  liquid  filled  projectiles, 3, 4  an<j 
eigenfrequencies  of  the  perturbed  fluid  motion,5’®  both  of  which  are  of  critical 
importance  in  determing  the  projectile's  stability. 
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t.  fc'fc  R 


NOMENCLATURE 


F  Froude  number,  (siR)Z  / $  H 

g  axial  acceleration  (  =  gravity  in  a  laboratory  frame) 

H  height  of  cylinder 

i  time  index,  tk  =  i  A  t 

j  radial  grid  index,  rj  =  (j-1)  Ar 

k  axial  grid  index,  zk  =  (k-1)  Az 

L  initial  fluid  level  before  spin-up 

Lz  dimensionless  fluid  angular  momentum  about  the  cylinder’s  axis,  Eq.  (A-9) 

NJ  number  of  radial  grid  points 

NK  number  of  axial  grid  points  used  in  Stage  3  free  surface  integrals 

p  pressure 

(r,  ©,  z)  dimensional  cylindrical  coordinates 

R  radius  of  cylinder 

r  r/R 

ro*  rH  dimensionless  radii  at  which  endwalls  are  intersected  (Figure  1) 

®FS  radial  coordinate  description  of  the  free  surface  during  Stage  3 

Re  Reynolds  number,  il  R2/V 

t  jxT 

ts  dimensionless  spin-up  time 

(u,v,w)  dimensionless  velocities  in  the  (r,  ©,  z)  directions,  respectively,  normalized 
by  -QR 

utel»  uBEL  contributions  to  u  from  the  top  and  bottom  Ekman  layers,  respectively 

UFS  contribution  to  u  from  free  surface  motion 

z  z/H 

Zps  axial  coordinate  description  of  the  free  surface  during  Stages  l  and  2b 

Z<>»  ^1  dimensionless  fluid  levels  at  the  axis  and  sidewall,  respectively 

relaxation  factor  in  velocity  iterations 
relaxation  factor  in  free  surface  iterations 
Ekman  boundary  layer  thickness 
kinematic  viscosity 
density 

angular  velocity  of  cylinder,  rad./sec. 
dimensionless  local  angular  velocity,  yr/r  *  <*>/■&. 

(  )  indicates  a  dimensional  variable 
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Appendix  A 

FREE  SURFACE  EQUATIONS 


In  this  appendix  we  develop  the  equations  used  to  describe  the  free  surface 
motion.  Following  Winch2  and  Goller  and  Ranov,2^  we  assume  that  surface  tension  can 
be  neglected,  in  which  case  balancing  the  centrifugal  and  gravitational  body  forces 
acting  on  a  fluid  element  at  the  surface  yields: 


cii 


-  2  - 
oj  r 


oL  r 


? 


% 


(A-l) 


or  in  dimensionless  form  (see  Nomenclature) 

ot2pS 
cLr 


F  uy  r  - 


F  r" 


(A-2) 


where  z  =  Zpg(r,  t)  defines  the  free  surface  contour.  This  equation  is  developed  and 
used  differently  depending  on  whether  the  flow  is  in  Stage  1,  2  or  3  (in  Gerber's22 
notation)  and  each  of  these  is  discussed  separately  below. 


Stage  1 

This  is  the  case  treated  by  Goller  and  Ranov,  in  which  neither  endwall  is 
intersected.  Let  Zj  =  Zps(l,  t),  i.e.,  the  dimensionless  fluid  level  at  the  sidewall. 
Then  from  (A-2), 

r1  z 

^rs(r,t)  -  zit)-Fj  —  dr'  (a-3) 

r  r 

Conservation  of  fluid  volume  requires  that 

znf  2PS  rdr  e  7T  ~  (A-4) 

o 

Substituting  (A-3)  into  (A-4)  gives 

ij'rdr  -  F^rdr/'  dr’}  -  7r  A 
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The  double  integration  in  the  second  term  can  be  simplified  through  an  integration  by 
parts  with  respect  to  the  dummy  variable  r'.*  The  result,  after  rearranging  and 
cancelling  common  factors,  is  the  following  expression  for  Z\i 


Z 


_L 

H 


+  F 


S 


v  r dr 


(A -5) 


Eqs.  (A-3)  and  (A-5)  together  define  the  free  surface  contour,  given  v(r,  t).  For 
completeness,  we  note  that  the  on-axis  fluid  level,  Z0  =  Zpg(0,  t)  can  be  inferred  to  be: 


TT-Ff'»‘(T  ~  r ) 


(A-6) 


If  the  integration  by  parts  noted  above  is  used  to  simplify  Goller  and  Ranov's  Eq.  (lc), 
the  result  is  consistent  with  (A-6). 


In  the  steady  limit  of  solid-body  rotation,  v  =  r,  and  the  integrations  in  Eqs.  (A-5) 
and  (A-6)  can  be  done  analytically.  One  can  easily  show  that 

L  _  i,(t  t  Z,(t  »«)  ,A_„X 

T  7  (A  7) 

From  this  one  can  infer  under  what  conditions  the  flow  will  remain  in  Stage  1.  For 
the  free  surface  to  first  intersect  the  bottom  endwall  (Stage  2a)  requires  that  Z0  =  0 
and  Zj  <  1,  for  which  (A-6)  and  (A-7)  together  imply  L/H  <  1/2  and  Fi  4  L/H.  To 
intersect  the  top  endwall  first  (Stage  2b)  requires  Z0  >  0  and  Z1  =  l,  for  which  (A-5) 
and  (A-7)  together  imply  L/H  >  1/2  and  F  >  4  (1  -  L/H).  For  present  purposes  we 
are  primarily  concerned  with  fill  ratios  greater  than  one-ha  ,  and  so  the  flow  will 
remain  in  Stage  1  for  all  time  if 

F  <  4  (l  -  A )  (A-8) 


As  a  reference  for  judging  the  degree  to  which  the  fluid  has  spun  up,  we  define 
a  dimensionless  angular  momentum  about  the  axis  as 


Lg(t) 


ZTtp  H  a  R4 


/' 


ZP$  vr* d  r 


(A-9) 


The  author  is  indebted  to  Mr.  Nathan  Gerber  for  pointing  this  out. 
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For  flows  remaining  in  Stage  1,  in  the  limit  of  solid-body  rotation,  Eqs.  (A-3), 
(A-5),  and  (A-9)  then  reduce  to 


L  (t  -  eo)  —(  —  *■  — 
r  A  v  H  1 2  ) 


(A-10) 


Stage  2b 


We  denote  by  rjj(t)  the  dimensionless  radius  at  which  the  top  endwall  is 
intersected.  To  the  right  of  this  point,  Zpg  =  1;  for  smaller  radii,  we  get  immediately 
from  Eq.  (A-2), 


t;  1  -  F 


Jr»  X* 


d  r' 


(A-ll) 


ZaW> 

o 


/*  rH  tr2 

,  -f  J  _  dr 


(A-12) 


Equation  (A-ll)  alone  is  insufficient  to  define  Zpg  given  v(r,  t),  since  rjj  is  also 
unknown.  An  additional  relation  for  rjj  can  be  derived  as  follows.  Replace  Eq.  (A-ll) 
with 


Z«<v,  t )  -  2 


(A— 13) 


and  again  invoke  conservation  of  fluid  volume: 


nf  "i^rdr  ♦  «  (,-rJ)  -  7T  L 


(A-14) 


Substituting  from  (A-13)  one  obtains 


27T  {*•/  "rdr  +l  "rdrl  •  nH 

Integrating  by  parts  again  simplifies  the  term  containing  the  double  integral,  and  the 
result  is  an  alternate  expression  for  Zn, 
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Equating  this  to  (A-12),  several  terms  cancel,  and  we  are  left  with 


As  shown  in  Section  2,  this  result  can  be  used  to  infer  ry,  and  then  Eq.  (A-ll)  used  for 
the  rest  of  the  contour. 

The  steady  limit  of  Eq.  (A-16)  yields 

(t  — )  -  [y  (>-  y)] 

and  Eq.  (A-12)  then  becomes 

[rC’-tt)]  '* 

Since  ZQ  >  0  for  the  flow  to  remain  in  Stage  2b,  we  see  from  Eqs.  (A-8)  and  (A-18)  that 
F  must  then  satisfy 

•»(>-  7r)  *  F  ‘  (’-  -w)"  (A-19) 

For  flows  in  this  stage,  the  solid-body  rotation  limit  of  Eq.  (A-9)  is  found  to  be 

-  T  (A'20) 

Though  not  treated  here,  results  analogous  to  these  can  also  be  derived  for  Stage 
2a;  the  steady-state  results  can  be  found  in  Reference  22. 

Stage  3 

For  the  reasons  discussed  in  the  main  text,  during  this  state  we  choose  z  as  our 
independent  variable  and  let  r  --  Rpg(z,  t)  describe  the  free  surface  contour.  Equation 
(A-2)  is  rewritten  as 


(A-17) 


(A-18) 
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where  now  it  is  implicitly  understood  that  to(z,  t)  =  o>(Rpg(z,  t),  t).  Defining  r0(t) 
=  Rfs(°>  0  and  integrating  (A-21)  from  the  bottom  up  gives 

*"“■*>  ■  r'<t>+  t/'z&FTo 

“  *  (A-22) 


-  r‘(t>  +  y  f' 

o 

■  **'  +  $/ 


RpS  (  z\t  )  oLz 
t) 


Conservation  of  mass  requires 


d  * 


Substituting  (A-22)  into  (A-23), 


-it 


cLi  ' 

co*(z' *) 


(A-23) 


The  second  term  can  again  be  reduced  to  a  single  integral  through  an  integration  by 
parts.  Resetting  ui  =  v/r,  the  result  is 

z...  L  2  r1  *!*<■* *tUf-M)de 


r  it)  *  i 
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_  _L  -  L  f 

H  F  *4  2rz  (  g ,  t  ) 


(A-24) 


where  of  course  v(z,  t)  is  understood  to  mean  v(Rps(z>  *)»  t).  Substitution  of  (A-24)  in 
(A-22)  yields 


r*  <f>  *  ’-T+jI 


i  n>  * 


z  dz 


(A-25) 


"  "  F  J0  m.t) 

In  principle,  Eqs.  (A-22)  and  (A-24)  define  the  free  surface  contour  for  a  given 

v  distribution.  But  Eq.  (A-22)  really  represents  an  integral  equation  for  the  desired 
2 

Rps>  and  one  which  cannot  be  inverted  in  closed  form.  As  discussed  in  Section  2,  an 
iterative  relaxation  procedure  was  used  to  solve  it. 

The  steady  limits  of  Eqs.  (A-24)  and  (A-25)  are 


,  r  l 

1  “l  ,/z 

mlm[—T 

“  F  ] 

(A-26) 

.  r  l 

i  1  ,/a 

*')  *[’-  H 

*  7J 

(A-27) 
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From  Eq.  (A-19)  we  infer  that  the  flow  will  enter  Stage  3  only  if 

f  ^  (t  -  jr)'1  (A_28) 

which,  through  (A-26),  insures  that  r0  >  0,  and,  through  (A-27),  that  r^  £  1.  Note  that 
in  the  limit  the  free  surface  forms  an  inner  cylinder  of  constant  radius,  rQ  =  rjj  = 

(1  -  L/H)l/2. 

Finally,  the  steady-state  limit  of  the  angular  momentum,  Eq.  (A-9),  in  this 
configuration  can  be  shown  to  be 

’  il’-O-jr)*- 7i*]  »-»> 

We  close  by  noting  that  all  the  expressions  given  above  for  the  limiting  condition 
of  solid-body  rotation  agree  with  those  derived  by  Gerber.22*23 
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Appendix  B 

FINITE  DIFFERENCE  EQUATIONS  FOR  THE  AZIMUTHAL  VELOCITY 

When  the  differencing  scheme  described  in  Section  2  is  applied  to  Eq.  (3),  the 
following  tridiagonal  system  of  equations  results: 


L *  1  1*1  1*1 

Ri  V  +  Bi  vi  +  c>  vi-' 
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Such  systems  are  efficiently  solved  via  a  special  adaptation  of  Gaussian  elimination 
(Reference  25).  Let 

vt’  '  £i  vl:’  *  Fi  ®-2) 


which,  when  substituted  into  Eq.  (B>1),  requires  that 

F): 

E  •  —  -  - ^  : - 

j  +  Cj  £,•-•> 


- 


D:  -  C :  F4-„, 


(B-3a) 


(B-3b) 
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These  are  applied  at  j  =  2,3  .  .  .  NJ  -  1,  with  the  starting  values 


E  =  F  -  o 

in  accordance  with  Eq.  (4a).  Then,  beginning  with  Eq.  (4b),  Eq.  (B-2)  is  used  to  obtain 
the  remaining  velocities. 

In  Stage  3,  the  above  procedure  is  modified  since  Eq.  (4a)  no  longer  applies.  Let 
rj*  denote  the  first  radial  grid  station  that  exceeds  r0;  then  only  the  points  j*  .  .  .  NJ 
will  be  included  in  the  algorithm  described  by  Eq.  (B-2)  and  (B-3).  The  starting  values 
for  E  and  F  are  inferred  by  requiring  that  a  linear  extrapolation  from  the  velocities  at 
rj*  and  rj*+i  to  rQ  should  satisfy  Eq.  (20).  The  results  are 

Ei*  *  (rik  ~r*)/(Ar  +  irj*  '  ) 

Fj*  ~  r  Ar y  (Ar  +  -  ra  ) 

The  values  used  for  the  relaxation  parameter  oc.  in  Eq.  (5)  are  given  in  the  table 

below. 


Table  B-l 

Values  of  Velocity  Relaxation  Parameter 


Free  Surface  Configuration  (Figure  1)  oo 


Z0  >  0.2 

1.0 

Stages  1,  2b 

0.1  1  Z0  <  0.2 

0.5 

Z0  <  0.1 

0.2 

Stage  3 

*0  >  0 

0.4 

These  values,  determined  from  experience,  are  not  necessarily  optimal.  But  with  them, 
convergence  was  achieved  typically  within  2-3  iterations. 
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